Application of Spatial Point Patterns Analysis to discover the geographical distribution of Grab hailing services in Singapore

Take-home_Ex01

Published

January 15, 2024

Modified

February 5, 2024

Setting the Scene

Human mobility, the movement of human beings in space and time, reflects the spatial-temporal characteristics of human behavior. With the advancement of Information and Communication Technologies (ICT) especially smart phone, a large volume of data related to human mobility have been collected. By using appropriate GIS analysis methods, these data are potentially useful in supporting smart city planning and management.

In Singapore, one of the important source of data related to human mobility is from Land Transport Authority (LTA) DataMall. Two data sets related to human mobility are provided by the portal, they are: Passenger Volume by Origin Destination Train Stations and Passenger Volume by Origin Destination Bus Stops. One of the limitation of these data sets is that their location are biased to either bus stops or MRT/LRT stations. In 2020, another very interesting human mobility data set called Grab Posisi was released by GRAB, one of the largest shared taxi operator in South-east Asia. There are two data sets been released and one of them is for Singapore.

Objectives

Geospatial analytics hold tremendous potential to address complex problems facing society. In this study, you are tasked to apply appropriate spatial point patterns analysis methods to discover the geographical and spatio-temporal distribution of Grab hailing services locations in Singapore.

The Task

  • Using appropriate function of sf and tidyverse, preparing the following geospatial data layer in sf tibble data.frames:

    • Grab taxi location points either by origins or destinations.

    • Road layer within Singapore excluding outer islands.

    • Singapore boundary layer excluding outer islands.

  • Using the extracted data, derive traditional Kernel Density Estimation layers

  • Using the extracted data, derive either Network Kernel Density Estimation (NKDE) or Temporal Network Kernel Density Estimation (TNKD)

  • Using appropriate tmap function, display the kernel density layers on OpenStreetMap of Singapore

  • Describe the spatial patterns revealed by the kernel density maps.

The Data

Aspatial data

For the purpose of this assignment, Grab-Posisi of Singapore will be used

Geospatial data

  • Road data set from OpenStreetMap of Geofabrik download server. The Malaysia, Singapore, and Brunei coverage should be downloaded.

  • Master Plan 2019 Subzone Boundary (No Sea) from Data.gov.sg

Installing and Loading the R packages

pacman::p_load(maptools, sf, raster, spatstat, tmap, tidyverse, arrow, spNetwork, classInt, viridis)

Data Wrangling

Importing Grab-Posisi Dataset

Since the dataset is stored as .parquet file format, we need to use the arrow package

grab_data <- arrow::open_dataset("../../data/aspatial/GrabPosisi") %>%
  collect()

Let’s convert the pingtimestamp column from integer to a datetime format

grab_data$pingtimestamp <- as_datetime(grab_data$pingtimestamp)

Extracting trip starting locations

origin_df <- grab_data %>%
  group_by(trj_id) %>%
  arrange(pingtimestamp) %>% # from dplyr package
  filter(row_number()==1) %>%
  mutate(weekday = wday(pingtimestamp, # functions from lubridate package
                        label=TRUE,
                        abbr=TRUE),
         start_hr = factor(hour(pingtimestamp)), # With factor(), it makes it ordinal
         day = factor(mday(pingtimestamp)))
view(origin_df)

Extracting trip destination locations

destination_df <- grab_data %>%
  group_by(trj_id) %>%
  arrange(desc(pingtimestamp)) %>% # from dplyr package
  filter(row_number()==1) %>%
  mutate(weekday = wday(pingtimestamp, # functions from lubridate package
                        label=TRUE,
                        abbr=TRUE),
         end_hr = factor(hour(pingtimestamp)), # With factor(), it makes it ordinal
         day = factor(mday(pingtimestamp)))
view(destination_df)

Converting into sf tibble dataframe

Let’s convert it into a sf tibble dataframe. According to the data source, our data is in Geodetic CRS: WGS84, so we need to convert it to SVY21 too.

origin_sf <- st_as_sf(origin_df, coords = c("rawlng", "rawlat"), crs = 4326) %>%
  st_transform(crs=3414)

destination_sf <- st_as_sf(destination_df, coords = c("rawlng", "rawlat"), crs = 4326) %>%
  st_transform(crs=3414)
st_geometry(origin_sf)
st_geometry(destination_sf)
head(origin_sf)
Simple feature collection with 6 features and 10 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 17201.85 ymin: 29236.55 xmax: 31518.05 ymax: 47967.97
Projected CRS: SVY21 / Singapore TM
# A tibble: 6 × 11
# Groups:   trj_id [6]
  trj_id driving_mode osname pingtimestamp        speed bearing accuracy weekday
  <chr>  <chr>        <chr>  <dttm>               <dbl>   <int>    <dbl> <ord>  
1 70895  car          andro… 2019-04-08 00:09:26  9.95      111      4   Mon    
2 21926  car          andro… 2019-04-08 00:09:48 11.0        75      4   Mon    
3 47498  car          ios    2019-04-08 00:09:50 18.6       307      8   Mon    
4 18103  car          andro… 2019-04-08 00:09:55  0.404     159      3   Mon    
5 41322  car          andro… 2019-04-08 00:09:57 17.9       232      3.9 Mon    
6 64813  car          ios    2019-04-08 00:10:03 17.1       106     10   Mon    
# ℹ 3 more variables: start_hr <fct>, day <fct>, geometry <POINT [m]>

Importing the Boundary of Singapore

mpsz_sf <- st_read(dsn ="../../data/geospatial/MP14_SUBZONE_WEB_PL", layer = "MP14_SUBZONE_WEB_PL") %>%
  st_transform(crs = 3414)
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\russelltjh\IS415-GAA\data\geospatial\MP14_SUBZONE_WEB_PL' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
sg_sf <- st_union(mpsz_sf)

plot(sg_sf)

st_geometry(sg_sf)
Geometry set for 1 feature 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21 / Singapore TM

Extracting Study Areas from Singapore Boundary

bedok <- mpsz_sf %>%
  filter(PLN_AREA_N == "BEDOK")
tamp <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
clementi <- mpsz_sf %>%
  filter(PLN_AREA_N == "CLEMENTI")
jurongwest <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")

Importing Road data from OpenStreetMap

We have the ESRI files downloaded, so let’s check what layers we can choose from this dataset.

file_path = "../../data/geospatial/malaysia-singapore-brunei.shp"
layers <- st_layers(file_path)
print(layers)
Driver: ESRI Shapefile 
Available layers:
                   layer_name geometry_type features fields crs_name
1  gis_osm_buildings_a_free_1       Polygon  1275019      5   WGS 84
2    gis_osm_landuse_a_free_1       Polygon    61509      4   WGS 84
3    gis_osm_natural_a_free_1       Polygon      665      4   WGS 84
4      gis_osm_natural_free_1         Point    45725      4   WGS 84
5     gis_osm_places_a_free_1       Polygon     1081      5   WGS 84
6       gis_osm_places_free_1         Point     6920      5   WGS 84
7       gis_osm_pofw_a_free_1       Polygon     5992      4   WGS 84
8         gis_osm_pofw_free_1         Point     2369      4   WGS 84
9       gis_osm_pois_a_free_1       Polygon    54173      4   WGS 84
10        gis_osm_pois_free_1         Point    80301      4   WGS 84
11    gis_osm_railways_free_1   Line String     5252      7   WGS 84
12       gis_osm_roads_free_1   Line String  1766163     10   WGS 84
13   gis_osm_traffic_a_free_1       Polygon    17502      4   WGS 84
14     gis_osm_traffic_free_1         Point    74074      4   WGS 84
15 gis_osm_transport_a_free_1       Polygon     1320      4   WGS 84
16   gis_osm_transport_free_1         Point    16991      4   WGS 84
17   gis_osm_waterways_free_1   Line String    34559      5   WGS 84
18     gis_osm_water_a_free_1       Polygon    19682      4   WGS 84

We only want the data for roads, so let’s pick the gis_osm_roads_free_1 layer.

MySgBr_Roads <- st_read(dsn ="../../data/geospatial/malaysia-singapore-brunei.shp", layer = "gis_osm_roads_free_1")

We notice that the MySgBr_Roads data is in Geodetic CRS: WGS 84 format, let’s convert it to SVY21.

MySgBr_Roads_SVY21 <- st_transform(MySgBr_Roads, crs = 3414)

Now that the data is in SVY21 format, let’s extract only the roads in Singapore using the Singapore boundary data.

indices <- st_intersects(sg_sf, MySgBr_Roads_SVY21)
sg_roads <- MySgBr_Roads_SVY21[unlist(indices), ]
print(sg_roads)
Simple feature collection with 233450 features and 10 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 2676.964 ymin: 15901.59 xmax: 56061.93 ymax: 50239.93
Projected CRS: SVY21 / Singapore TM
First 10 features:
     osm_id code      fclass              name  ref oneway maxspeed layer
1   4386520 5113     primary      Orchard Road <NA>      F       50     0
18  4887867 5122 residential  Hougang Avenue 1 <NA>      B       50     0
497 8096835 5113     primary       Scotts Road <NA>      F       60     0
505 9584642 5115    tertiary     Keng Lee Road <NA>      F       50     0
506 9584847 5153     footway              <NA> <NA>      B        0     0
507 9585045 5113     primary       Newton Road <NA>      F       60     0
508 9585074 5122 residential      Sarkies Road <NA>      B       50     0
509 9585621 5113     primary     Paterson Road <NA>      F       50     0
510 9585771 5113     primary Orchard Boulevard <NA>      F       50     0
511 9586040 5113     primary     Paterson Road <NA>      F       50     0
    bridge tunnel                       geometry
1        F      F LINESTRING (27637.52 32038....
18       F      F LINESTRING (34016.12 36785....
497      F      F LINESTRING (28588.71 32765....
505      F      F LINESTRING (29165.25 32888....
506      F      F LINESTRING (28834.55 32990....
507      F      F LINESTRING (28664.35 32858....
508      F      F LINESTRING (28433.13 32988....
509      F      F LINESTRING (27835.44 31930....
510      F      F LINESTRING (28168.89 31418....
511      F      F LINESTRING (27685.73 31733....

Filtering Road types

Let’s check the road classes in the dataset.

# Assuming sg_roads is your dataset
unique_road_types <- sg_roads %>% 
  distinct(fclass) %>% 
  arrange(fclass) # Optional: to sort the road types alphabetically

# Print the unique road types
print(unique_road_types)
           fclass
1       bridleway
2        cycleway
3         footway
4   living_street
5        motorway
6   motorway_link
7            path
8      pedestrian
9         primary
10   primary_link
11    residential
12      secondary
13 secondary_link
14        service
15          steps
16       tertiary
17  tertiary_link
18          track
19   track_grade1
20   track_grade2
21   track_grade3
22   track_grade4
23   track_grade5
24          trunk
25     trunk_link
26   unclassified
27        unknown

We can see that the classes contains road types that does not allows car to be driven on, such as footway, cycleway, pedestrian and etc. We will make some assumptions on which classes do not allow cars and filter them out.

# Filter out road types not meant for cars
sg_car_roads <- sg_roads %>% 
  filter(!fclass %in% c("bridleway", "cycleway", "footway", "path", "pedestrian", 
                        "steps", "track", "track_grade1", "track_grade2", 
                        "track_grade3", "track_grade4", "track_grade5")) 

# View the filtered dataset
print(sg_car_roads)
Simple feature collection with 164952 features and 10 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 2676.964 ymin: 19123.75 xmax: 55604.55 ymax: 50239.93
Projected CRS: SVY21 / Singapore TM
First 10 features:
    osm_id code      fclass              name  ref oneway maxspeed layer bridge
1  4386520 5113     primary      Orchard Road <NA>      F       50     0      F
2  4887867 5122 residential  Hougang Avenue 1 <NA>      B       50     0      F
3  8096835 5113     primary       Scotts Road <NA>      F       60     0      F
4  9584642 5115    tertiary     Keng Lee Road <NA>      F       50     0      F
5  9585045 5113     primary       Newton Road <NA>      F       60     0      F
6  9585074 5122 residential      Sarkies Road <NA>      B       50     0      F
7  9585621 5113     primary     Paterson Road <NA>      F       50     0      F
8  9585771 5113     primary Orchard Boulevard <NA>      F       50     0      F
9  9586040 5113     primary     Paterson Road <NA>      F       50     0      F
10 9587072 5113     primary River Valley Road <NA>      F       50     0      F
   tunnel                       geometry
1       F LINESTRING (27637.52 32038....
2       F LINESTRING (34016.12 36785....
3       F LINESTRING (28588.71 32765....
4       F LINESTRING (29165.25 32888....
5       F LINESTRING (28664.35 32858....
6       F LINESTRING (28433.13 32988....
7       F LINESTRING (27835.44 31930....
8       F LINESTRING (28168.89 31418....
9       F LINESTRING (27685.73 31733....
10      F LINESTRING (28893.45 30872....

Extracting the Roads of Study Areas

indices <- st_intersects(bedok, sg_car_roads)
bedok_roads <- sg_car_roads[unlist(indices), ]

indices <- st_intersects(tamp, sg_car_roads)
tamp_roads <- sg_car_roads[unlist(indices), ]

indices <- st_intersects(clementi, sg_car_roads)
clementi_roads <- sg_car_roads[unlist(indices), ]

indices <- st_intersects(jurongwest, sg_car_roads)
jurongwest_roads <- sg_car_roads[unlist(indices), ]
print(bedok_roads)
Simple feature collection with 8080 features and 10 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 32694.22 ymin: 31745.24 xmax: 42410.18 ymax: 36884.38
Projected CRS: SVY21 / Singapore TM
First 10 features:
       osm_id code        fclass                 name  ref oneway maxspeed
506  22698750 5113       primary   Marine Parade Road <NA>      F       60
514  22699682 5113       primary Bedok South Avenue 1 <NA>      F       60
516  22699687 5131 motorway_link BEDOK SOUTH AVENUE 1 <NA>      F       50
517  22699689 5131 motorway_link                 <NA> <NA>      F       50
518  22699690 5131 motorway_link Bedok South Avenue 1 <NA>      F       50
519  22699691 5131 motorway_link Bedok South Avenue 1 <NA>      F       50
1668 23825774 5131 motorway_link                 <NA> <NA>      F       50
1799 23924790 5141       service                 <NA> <NA>      F        0
1818 23960766 5115      tertiary         Marine Vista <NA>      B       50
1819 23960768 5131 motorway_link                 <NA> <NA>      F       50
     layer bridge tunnel                       geometry
506      0      F      F LINESTRING (39026.39 32545....
514      0      F      F LINESTRING (39247.54 32654....
516      1      T      F LINESTRING (39326.01 32324....
517      0      F      F LINESTRING (39346.99 32375....
518      1      T      F LINESTRING (39318.08 32396....
519      0      F      F LINESTRING (39336.16 32327....
1668     0      F      F LINESTRING (38343.1 32057.5...
1799     0      F      F LINESTRING (37847.67 32166....
1818     0      F      F LINESTRING (37771.93 31948....
1819     0      F      F LINESTRING (37771.93 31948....

Displaying Road Study Area

ggplot(data = clementi_roads) +
  geom_sf() +
  theme_minimal() +
  labs(title = "Roads in Clementi", 
       x = "Longitude", 
       y = "Latitude")

Data Preparation

Creating ppp object

spatstat requires the analytical data in ppp object form. Let’s convert our grab_sf into a ppp object.

origin_ppp <- as.ppp(origin_sf)
destination_ppp <- as.ppp(destination_sf)

Now, we will display the properties of our new ppp object.

summary(origin_ppp)
Marked planar point pattern:  28000 points
Average intensity 2.473666e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

marks are of type 'character'
Summary:
   Length     Class      Mode 
    28000 character character 

Window: rectangle = [3628.24, 49845.23] x [25198.14, 49689.64] units
                    (46220 x 24490 units)
Window area = 1131920000 square units
summary(destination_ppp)
Marked planar point pattern:  28000 points
Average intensity 2.493661e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

marks are of type 'character'
Summary:
   Length     Class      Mode 
    28000 character character 

Window: rectangle = [3637.21, 49870.63] x [25221.3, 49507.79] units
                    (46230 x 24290 units)
Window area = 1122850000 square units

Creating owin object

It is good practice to confine the analysis within a geographical area like Singapore’s boundaries. In spatstat, we can use the owin object to represent this polygonal region.

sg_owin <- as.owin(sg_sf)
bedok_owin <- as.owin(bedok)
tamp_owin <- as.owin(tamp)
clementi_owin <- as.owin(clementi)
jurongwest_owin <- as.owin(jurongwest)

Combining point events object and owin object

originSG_ppp = origin_ppp[sg_owin]
originBedok_ppp = origin_ppp[bedok_owin]
originTamp_ppp = origin_ppp[tamp_owin]
originClementi_ppp = origin_ppp[clementi_owin]
originJurongWest_ppp = origin_ppp[jurongwest_owin]

destinationSG_ppp = destination_ppp[sg_owin]
destinationBedok_ppp = destination_ppp[bedok_owin]
destinationTamp_ppp = destination_ppp[tamp_owin]
destinationClementi_ppp = destination_ppp[clementi_owin]
destinationJurongWest_ppp = destination_ppp[jurongwest_owin]
summary(originSG_ppp)
Marked planar point pattern:  28000 points
Average intensity 3.580815e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

marks are of type 'character'
Summary:
   Length     Class      Mode 
    28000 character character 

Window: polygonal boundary
80 separate polygons (35 holes)
                  vertices         area relative.area
polygon 1            14650  6.97996e+08      8.93e-01
polygon 2 (hole)         3 -2.21090e+00     -2.83e-09
polygon 3              285  1.61128e+06      2.06e-03
polygon 4 (hole)         3 -2.05920e-03     -2.63e-12
polygon 5 (hole)         3 -8.83647e-03     -1.13e-11
polygon 6              668  5.40368e+07      6.91e-02
polygon 7               44  2.26577e+03      2.90e-06
polygon 8               27  1.50315e+04      1.92e-05
polygon 9              711  1.28815e+07      1.65e-02
polygon 10 (hole)       36 -4.01660e+04     -5.14e-05
polygon 11 (hole)      317 -5.11280e+04     -6.54e-05
polygon 12 (hole)        3 -3.41405e-01     -4.37e-10
polygon 13 (hole)        3 -2.89050e-05     -3.70e-14
polygon 14              77  3.29939e+05      4.22e-04
polygon 15              30  2.80002e+04      3.58e-05
polygon 16 (hole)        3 -2.83151e-01     -3.62e-10
polygon 17              71  8.18750e+03      1.05e-05
polygon 18 (hole)        3 -1.68316e-04     -2.15e-13
polygon 19 (hole)       36 -7.79904e+03     -9.97e-06
polygon 20 (hole)        4 -2.05611e-02     -2.63e-11
polygon 21 (hole)        3 -2.18000e-06     -2.79e-15
polygon 22 (hole)        3 -3.65501e-03     -4.67e-12
polygon 23 (hole)        3 -4.95057e-02     -6.33e-11
polygon 24 (hole)        3 -3.99521e-02     -5.11e-11
polygon 25 (hole)        3 -6.62377e-01     -8.47e-10
polygon 26 (hole)        3 -2.09065e-03     -2.67e-12
polygon 27              91  1.49663e+04      1.91e-05
polygon 28 (hole)       26 -1.25665e+03     -1.61e-06
polygon 29 (hole)      349 -1.21433e+03     -1.55e-06
polygon 30 (hole)       20 -4.39069e+00     -5.62e-09
polygon 31 (hole)       48 -1.38338e+02     -1.77e-07
polygon 32 (hole)       28 -1.99862e+01     -2.56e-08
polygon 33              40  1.38607e+04      1.77e-05
polygon 34 (hole)       40 -6.00381e+03     -7.68e-06
polygon 35 (hole)        7 -1.40545e-01     -1.80e-10
polygon 36 (hole)       12 -8.36709e+01     -1.07e-07
polygon 37              45  2.51218e+03      3.21e-06
polygon 38             142  3.22293e+03      4.12e-06
polygon 39             148  3.10395e+03      3.97e-06
polygon 40              75  1.73526e+04      2.22e-05
polygon 41              83  5.28920e+03      6.76e-06
polygon 42             211  4.70521e+05      6.02e-04
polygon 43             106  3.04104e+03      3.89e-06
polygon 44             266  1.50631e+06      1.93e-03
polygon 45              71  5.63061e+03      7.20e-06
polygon 46              10  1.99717e+02      2.55e-07
polygon 47             478  2.06120e+06      2.64e-03
polygon 48             155  2.67502e+05      3.42e-04
polygon 49            1027  1.27782e+06      1.63e-03
polygon 50 (hole)        3 -1.16959e-03     -1.50e-12
polygon 51              65  8.42861e+04      1.08e-04
polygon 52              47  3.82087e+04      4.89e-05
polygon 53               6  4.50259e+02      5.76e-07
polygon 54             132  9.53357e+04      1.22e-04
polygon 55 (hole)        3 -3.23310e-04     -4.13e-13
polygon 56               4  2.69313e+02      3.44e-07
polygon 57 (hole)        3 -1.46474e-03     -1.87e-12
polygon 58            1045  4.44510e+06      5.68e-03
polygon 59              22  6.74651e+03      8.63e-06
polygon 60              64  3.43149e+04      4.39e-05
polygon 61 (hole)        3 -1.98390e-03     -2.54e-12
polygon 62 (hole)        4 -1.13774e-02     -1.46e-11
polygon 63              14  5.86546e+03      7.50e-06
polygon 64              95  5.96187e+04      7.62e-05
polygon 65 (hole)        4 -1.86410e-02     -2.38e-11
polygon 66 (hole)        3 -5.12482e-03     -6.55e-12
polygon 67 (hole)        3 -1.96410e-03     -2.51e-12
polygon 68 (hole)        3 -5.55856e-03     -7.11e-12
polygon 69             234  2.08755e+06      2.67e-03
polygon 70              10  4.90942e+02      6.28e-07
polygon 71             234  4.72886e+05      6.05e-04
polygon 72 (hole)       13 -3.91907e+02     -5.01e-07
polygon 73              15  4.03300e+04      5.16e-05
polygon 74             227  1.10308e+06      1.41e-03
polygon 75              10  6.60195e+03      8.44e-06
polygon 76              19  3.09221e+04      3.95e-05
polygon 77             145  9.61782e+05      1.23e-03
polygon 78              30  4.28933e+03      5.49e-06
polygon 79              37  1.29481e+04      1.66e-05
polygon 80               4  9.47108e+01      1.21e-07
enclosing rectangle: [2667.54, 56396.44] x [15748.72, 50256.33] units
                     (53730 x 34510 units)
Window area = 781945000 square units
Fraction of frame area: 0.422
summary(destinationSG_ppp)
Marked planar point pattern:  27997 points
Average intensity 3.580432e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

marks are of type 'character'
Summary:
   Length     Class      Mode 
    27997 character character 

Window: polygonal boundary
80 separate polygons (35 holes)
                  vertices         area relative.area
polygon 1            14650  6.97996e+08      8.93e-01
polygon 2 (hole)         3 -2.21090e+00     -2.83e-09
polygon 3              285  1.61128e+06      2.06e-03
polygon 4 (hole)         3 -2.05920e-03     -2.63e-12
polygon 5 (hole)         3 -8.83647e-03     -1.13e-11
polygon 6              668  5.40368e+07      6.91e-02
polygon 7               44  2.26577e+03      2.90e-06
polygon 8               27  1.50315e+04      1.92e-05
polygon 9              711  1.28815e+07      1.65e-02
polygon 10 (hole)       36 -4.01660e+04     -5.14e-05
polygon 11 (hole)      317 -5.11280e+04     -6.54e-05
polygon 12 (hole)        3 -3.41405e-01     -4.37e-10
polygon 13 (hole)        3 -2.89050e-05     -3.70e-14
polygon 14              77  3.29939e+05      4.22e-04
polygon 15              30  2.80002e+04      3.58e-05
polygon 16 (hole)        3 -2.83151e-01     -3.62e-10
polygon 17              71  8.18750e+03      1.05e-05
polygon 18 (hole)        3 -1.68316e-04     -2.15e-13
polygon 19 (hole)       36 -7.79904e+03     -9.97e-06
polygon 20 (hole)        4 -2.05611e-02     -2.63e-11
polygon 21 (hole)        3 -2.18000e-06     -2.79e-15
polygon 22 (hole)        3 -3.65501e-03     -4.67e-12
polygon 23 (hole)        3 -4.95057e-02     -6.33e-11
polygon 24 (hole)        3 -3.99521e-02     -5.11e-11
polygon 25 (hole)        3 -6.62377e-01     -8.47e-10
polygon 26 (hole)        3 -2.09065e-03     -2.67e-12
polygon 27              91  1.49663e+04      1.91e-05
polygon 28 (hole)       26 -1.25665e+03     -1.61e-06
polygon 29 (hole)      349 -1.21433e+03     -1.55e-06
polygon 30 (hole)       20 -4.39069e+00     -5.62e-09
polygon 31 (hole)       48 -1.38338e+02     -1.77e-07
polygon 32 (hole)       28 -1.99862e+01     -2.56e-08
polygon 33              40  1.38607e+04      1.77e-05
polygon 34 (hole)       40 -6.00381e+03     -7.68e-06
polygon 35 (hole)        7 -1.40545e-01     -1.80e-10
polygon 36 (hole)       12 -8.36709e+01     -1.07e-07
polygon 37              45  2.51218e+03      3.21e-06
polygon 38             142  3.22293e+03      4.12e-06
polygon 39             148  3.10395e+03      3.97e-06
polygon 40              75  1.73526e+04      2.22e-05
polygon 41              83  5.28920e+03      6.76e-06
polygon 42             211  4.70521e+05      6.02e-04
polygon 43             106  3.04104e+03      3.89e-06
polygon 44             266  1.50631e+06      1.93e-03
polygon 45              71  5.63061e+03      7.20e-06
polygon 46              10  1.99717e+02      2.55e-07
polygon 47             478  2.06120e+06      2.64e-03
polygon 48             155  2.67502e+05      3.42e-04
polygon 49            1027  1.27782e+06      1.63e-03
polygon 50 (hole)        3 -1.16959e-03     -1.50e-12
polygon 51              65  8.42861e+04      1.08e-04
polygon 52              47  3.82087e+04      4.89e-05
polygon 53               6  4.50259e+02      5.76e-07
polygon 54             132  9.53357e+04      1.22e-04
polygon 55 (hole)        3 -3.23310e-04     -4.13e-13
polygon 56               4  2.69313e+02      3.44e-07
polygon 57 (hole)        3 -1.46474e-03     -1.87e-12
polygon 58            1045  4.44510e+06      5.68e-03
polygon 59              22  6.74651e+03      8.63e-06
polygon 60              64  3.43149e+04      4.39e-05
polygon 61 (hole)        3 -1.98390e-03     -2.54e-12
polygon 62 (hole)        4 -1.13774e-02     -1.46e-11
polygon 63              14  5.86546e+03      7.50e-06
polygon 64              95  5.96187e+04      7.62e-05
polygon 65 (hole)        4 -1.86410e-02     -2.38e-11
polygon 66 (hole)        3 -5.12482e-03     -6.55e-12
polygon 67 (hole)        3 -1.96410e-03     -2.51e-12
polygon 68 (hole)        3 -5.55856e-03     -7.11e-12
polygon 69             234  2.08755e+06      2.67e-03
polygon 70              10  4.90942e+02      6.28e-07
polygon 71             234  4.72886e+05      6.05e-04
polygon 72 (hole)       13 -3.91907e+02     -5.01e-07
polygon 73              15  4.03300e+04      5.16e-05
polygon 74             227  1.10308e+06      1.41e-03
polygon 75              10  6.60195e+03      8.44e-06
polygon 76              19  3.09221e+04      3.95e-05
polygon 77             145  9.61782e+05      1.23e-03
polygon 78              30  4.28933e+03      5.49e-06
polygon 79              37  1.29481e+04      1.66e-05
polygon 80               4  9.47108e+01      1.21e-07
enclosing rectangle: [2667.54, 56396.44] x [15748.72, 50256.33] units
                     (53730 x 34510 units)
Window area = 781945000 square units
Fraction of frame area: 0.422

Kernel Density Estimation

kde_originSG_bw <- density(originSG_ppp,
                           sigma=bw.diggle,
                           edge=TRUE,
                           kernel="gaussian")

plot(kde_originSG_bw)

bw <- bw.diggle(originSG_ppp)
bw
   sigma 
8.751942 
kde_destinationSG_bw <- density(destinationSG_ppp,
                                sigma=bw.diggle,
                                edge=TRUE,
                                kernel="gaussian")
plot(kde_destinationSG_bw)

bw <- bw.diggle(destinationSG_ppp)
bw
   sigma 
8.752411 

As seen above, there is a concentration of grab pick up and drop of at the airport area. We can assume that many travellers use grab much more than Singaporeans in their day to day life.

Let’s try using a fixed bandwidth to see if there’s any significant difference.

kde_destinationSG_bw <- density(destinationSG_ppp,
                                sigma=0.5,
                                edge=TRUE,
                                kernel="gaussian")
plot(kde_destinationSG_bw)

bw <- bw.diggle(destinationSG_ppp)
bw
   sigma 
8.752411 

Let’s try to see KDE in the Clementi Study Area

kde_originClementi_bw <- density(originClementi_ppp,
                           sigma=bw.diggle,
                           edge=TRUE,
                           kernel="gaussian")

plot(kde_originClementi_bw)

bw <- bw.diggle(originClementi_ppp)
bw
   sigma 
19.17829 
# Check the class of the KDE result
print(class(kde_originClementi_bw))
[1] "im"
# Print a summary of the KDE result
print(summary(kde_originClementi_bw))
real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [18445.49, 21581.52] x [30449.66, 35936.17] units
dimensions of each pixel: 24.5 x 42.86335 units
Image is defined on a subset of the rectangular grid
Subset area = 9517629.01527276 square units
Subset area fraction = 0.553
Pixel values (inside window):
    range = [-3.417405e-19, 0.003180279]
    integral = 494.2188
    mean = 5.192667e-05
# Convert 'im' object to a RasterLayer
kde_raster <- raster(kde_originClementi_bw)

# Set the CRS for the raster using EPSG code for SVY21/Singapore TM
crs(kde_raster) <- CRS("EPSG:3414")

# Create a base viridis palette
palette <- viridis(100, begin = 0, end = 1, direction = -1) # A smooth gradient
# Generating a palette with 100 colors, with Beginning, Ending and Direction

tmap_mode("view")

# Plotting with tmap using the revised approach
tm_shape(kde_raster) +
  tm_raster(alpha = 0.75, style = "cont", palette = palette, n = 5) +
  tm_basemap("OpenStreetMap") +
  tm_layout(main.title = "KDE Overlay on OpenStreetMap")
Customizing Viridis Palette

The customization involves specifying the number of colors (viridis(100, …)) and potentially adjusting the portion of the color map used (begin, end, and direction parameters). These adjustments allow for finer control over the visual representation of your data:

  • Number of Colors: By generating a palette with 100 colors, you create a more granular color transition in the overlay. This can help in distinguishing subtle differences in the KDE values across the map.

  • Begin, End, and Direction: These parameters offer control over which portion of the viridis color map is utilized and the direction of the color gradient. While not demonstrated explicitly in the simplified example, adjusting these parameters can tailor the visual impact of the overlay, emphasizing specific ranges of data values more effectively.

Re-scaling KDE Values

In the code below, we will re-scale the unit of measurement from meter to kilometer

originSG_ppp.km <- rescale(originSG_ppp, 1000, 'km')
destinationSG_ppp.km <- rescale(destinationSG_ppp, 1000, 'km')
kde_originSG.bw <- density(originSG_ppp.km,
                           sigma=bw.diggle,
                           edge=TRUE,
                           kernel="quartic")
plot(kde_originSG.bw)

kde_destinationSG.bw <- density(destinationSG_ppp.km,
                                sigma=bw.diggle,
                                edge=TRUE,
                                kernel="gaussian")
plot(kde_destinationSG.bw)

originBedok_ppp.km <- rescale(originBedok_ppp,1000,'km')
kde_originBedok.bw <- density(originBedok_ppp.km,
                              sigma=bw.diggle,
                              edge=TRUE,
                              kernel="quartic")
plot(kde_originBedok.bw)

Network Constrained Kernel Density Estimation (NKDE)

Clementi Area

First, let’s filter the origin and destination points that are within the Clementi area only.

# Filter origin points in Clementi
intersects_index <- st_intersects(origin_sf, clementi, sparse = FALSE)

# Assuming `intersects_index` is a logical matrix from st_intersects
intersects_vector <- apply(intersects_index, 1, any)

# Use this logical vector to subset origin_sf
origin_sf_filtered <- origin_sf[intersects_vector, ]
# Filter destination points in Clementi
intersects_index <- st_intersects(destination_sf, clementi, sparse = FALSE)

# Assuming `intersects_index` is a logical matrix from st_intersects
intersects_vector <- apply(intersects_index, 1, any)

# Use this logical vector to subset origin_sf
destination_sf_filtered <- destination_sf[intersects_vector, ]
# Lixelize the Roads Network
origin_lixels <- lixelize_lines(clementi_roads, 750, mindist = 375)

# Generate Line Center Points
origin_samples <- lines_center(origin_lixels)

# Conduct Network Constrained KDE 
origin_densities <- nkde(clementi_roads,
                  events = origin_sf_filtered,
                  w = rep(1, nrow(origin_sf_filtered)),
                  samples = origin_samples,
                  kernel_name = "quartic",
                  bw = 300, # Confirm this bandwidth is appropriate for your data
                  div = "bw",
                  method = "simple",
                  digits = 1,
                  tol = 1,
                  grid_shape = c(1, 1),
                  max_depth = 8,
                  agg = 5, # Aggregates events within a 5m radius
                  sparse = TRUE,
                  verbose = FALSE)

# Assign Densities and Rescale
origin_samples$density <- origin_densities * 1000
origin_lixels$density <- origin_densities * 1000

# Visualize with tmap
tmap_mode("view")
tm_shape(origin_lixels) +
  tm_lines(col = "density") +
tm_shape(origin_sf_filtered) +
  tm_dots()
tmap_mode('plot')
# Lixelize the Roads Network
destination_lixels <- lixelize_lines(clementi_roads, 750, mindist = 375)

# Generate Line Center Points
destination_samples <- lines_center(destination_lixels)

# Conduct Network Constrained KDE 
destination_densities <- nkde(clementi_roads,
                  events = destination_sf_filtered,
                  w = rep(1, nrow(destination_sf_filtered)),
                  samples = destination_samples,
                  kernel_name = "quartic",
                  bw = 300, # Confirm this bandwidth is appropriate for your data
                  div = "bw",
                  method = "simple",
                  digits = 1,
                  tol = 1,
                  grid_shape = c(1, 1),
                  max_depth = 8,
                  agg = 5, # Aggregates events within a 5m radius
                  sparse = TRUE,
                  verbose = FALSE)

# Assign Densities and Rescale
destination_samples$density <- destination_densities * 1000
destination_lixels$density <- destination_densities * 1000

# Visualize with tmap
tmap_mode("view")
tm_shape(destination_lixels) +
  tm_lines(col = "density") +
tm_shape(destination_sf_filtered) +
  tm_dots()
tmap_mode('plot')

Visualizing Origin

# Visualize with tmap
tmap_mode("view")
tm_shape(origin_lixels) +
  tm_lines(col = "density", scale = 2) + 
  tm_basemap("OpenStreetMap") +
  tm_layout(main.title = "NKDE: Grab Pick Up Locations in Clementi Area")
tmap_mode('plot')

The map above shows that the popular locations for grab pick ups are concentrated in roughly 3 areas:

  • “West Coast Road” along the HDB blocks, Wespeak Terrace, Westview Gardens, and other private housing nearby.

  • “Commonwealth Avenue West”, along Block 455 bus stop.

  • “Clementi Avenue 6”, along Clementi View, Clementi Ridges and Nan Hua Primary School.

The common trait that all these areas have is that they are residential areas.

Visualizing Destination

# Visualize with tmap
tmap_mode("view")
tm_shape(destination_lixels) +
  tm_lines(col = "density", scale = 2) + 
  tm_basemap("OpenStreetMap") +
  tm_layout(main.title = "NKDE: Grab Pick Up Locations in Clementi Area")
tmap_mode('plot')
Back to top